#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <mpi.h>
#include <gsl/gsl_math.h>

#include "halofinder.h"

float halo_truncation_radius(double m)
{
  double critdens;
  
  if(strcmp(haloFinderData.HaloTruncRadType,"MEANOVERDENS") == 0)
    return (float) (pow(m/(haloFinderData.HaloOverdens*haloFinderData.OmegaM*RHO_CRIT*4.0/3.0*M_PI),1.0/3.0));
  else if(strcmp(haloFinderData.HaloTruncRadType,"CRITOVERDENS") == 0)
    {
      critdens = (double) (RHO_CRIT*(haloFinderData.OmegaM/pow(haloFinderData.ScaleFactorSnapshot,3.0) 
				    + 1.0-haloFinderData.OmegaM)*pow(haloFinderData.ScaleFactorSnapshot,3.0));
      return (float) (pow(m/(haloFinderData.HaloOverdens*critdens*4.0/3.0*M_PI),1.0/3.0));
    }
  else
    {
      fprintf(stderr,"%d: No halo truncation radius rule or rule does not match one of options!\n",ThisTask);
      MPI_Abort(MPI_COMM_WORLD,123);
      return -1.0;
    }
}

float halo_exlcusion_radius(double m1, double m2)
{
  float r1,r2;
    
  if(strcmp(haloFinderData.HaloExclType,"HARDSPHERE") == 0)
    {
      r1 = halo_truncation_radius(m1);
      r2 = halo_truncation_radius(m2);
      return r1 + r2;
    }
  else if(strcmp(haloFinderData.HaloExclType,"SOFTSPHERE") == 0 || strcmp(haloFinderData.HaloExclType,"SO") == 0)
    {
      r1 = halo_truncation_radius(m1);
      r2 = halo_truncation_radius(m2);
      if(r1 > r2)
	return r1;
      else
	return r2;
    }
  else
    {
      fprintf(stderr,"%d: No halo exclusion rule or rule does not match one of options!\n",ThisTask);
      MPI_Abort(MPI_COMM_WORLD,123);
      return -1.0;
    }
}
